library(rgdal)
library(foreign)
library(data.table)
library(rhandsontable)
library(spdep)
library(tmap)
# Load districts and attributes
g2 <- readOGR("./data", "svyMaps_2016.06.22_sara")
OGR data source with driver: ESRI Shapefile 
Source: "./data", layer: "svyMaps_2016.06.22_sara"
with 2078 features
It has 11 fields
dt2 <- read.dta("./data/SSApoverty_Dist_forGWR.12.dta")
# Keep STATA labels
dt2.lbl <- data.table(varCode=names(dt2), varLabel=attr(dt2, "var.labels"))
setkey(dt2.lbl, varCode)
dt2 <- data.table(dt2)
# Merge attributes to shapefile
g2.dt <- data.table(g2@data)
g2.dt[, rn22 := as.character(rn)]
g2.dt[, rn := row.names(g2)]
# Recode Ethiopia woredas
dt2[svyCode=="eth2010", svyL2Cd := svyL1Cd * 10000 + svyL2Cd]
setkey(g2.dt, ISO3, svyCode, svyL1Cd, svyL2Cd)
setkey(dt2, ISO3, svyCode, svyL1Cd, svyL2Cd)
# Any unmatched?
dt2[!g2.dt][, .N, by=svyCode]
Empty data.table (0 rows) of 2 cols: svyCode,N
# Remove `rn` from Sara's vars
dt2[, rn := NULL]
# Seems okay, so merge
g2.dt <- dt2[g2.dt]
names(g2.dt)[names(g2.dt) %like% "rn"]
[1] "rn"   "rn22"
# Neighbor list
coords <- coordinates(g2)
nb2 <- poly2nb(g2, row.names=paste(g2$ISO3, g2$rn, sep="."))
nb2
Neighbour list object:
Number of regions: 2078 
Number of nonzero links: 10902 
Percentage nonzero weights: 0.2524731 
Average number of links: 5.246391 
3 regions with no links:
MOZ.424 MWI.757 MWI.780

Plot contiguities for e.g. Ghana

# Plot contiguities (e.g. GHA)
gha <- g2[g2$ISO3=="GHA",]
gha.dt <- g2.dt[ISO3=="GHA"]
gha <- SpatialPolygonsDataFrame(gha, data.frame(gha.dt), match.ID="rn")
coords.gha <- coordinates(gha)
nb2.gha <- poly2nb(gha)
plot(gha, border="grey",
  main=list("Contiguity - Ghana (170 districts)"))
plot(nb2.gha, coords.gha, col="blue", add=T)
legend("bottomleft", 
  legend=capture.output(summary(nb2.gha))[1:7], cex=.6)


# Spatial weights for SSA (and Ghana only)
#w2 <- nb2listw(nb2)
w.gha <- nb2listw(nb2.gha)

# Explanatory vars
var <- c("spei_lt", "L1_speihishock", "L1_speiloshock")

# Plot to identify districts with missing values (Ghana only)
tmap_mode("view")
qtm(gha, fill=var[1])

# Impute missing values with regional median
tmp <- gha.dt[, lapply(.SD, median, na.rm=T), by=svyL1Cd, .SDcols=var]
setnames(tmp, 2:4, paste(names(tmp), "imp", sep="_"))
setkey(tmp, svyL1Cd)
setkey(gha.dt, svyL1Cd)
gha.dt <- tmp[gha.dt]

gha.dt[is.na(spei_lt), spei_lt := spei_lt_imp]
gha.dt[is.na(L1_speihishock), L1_speihishock := L1_speihishock_imp]
gha.dt[is.na(L1_speiloshock), L1_speiloshock := L1_speiloshock_imp]
# Spatial weights for SSA (and Ghana only)
#w2 <- nb2listw(nb2)
w.gha <- nb2listw(nb2.gha)
# Explanatory vars
var <- c("spei_lt", "L1_speihishock", "L1_speiloshock")
# Plot to identify districts with missing values (Ghana only)
tmap_mode("view")
tmap mode set to interactive viewing
qtm(gha, fill=var[1], colorNA="grey", textNA="Missing")
# Impute missing values with regional median
tmp <- gha.dt[, lapply(.SD, median, na.rm=T), by=svyL1Cd, .SDcols=var]
setnames(tmp, paste(names(tmp), "imp", sep="_"))
setkey(tmp, svyL1Cd_imp)
setkey(gha.dt, svyL1Cd)
gha.dt <- tmp[gha.dt]
gha.dt[is.na(spei_lt), spei_lt := spei_lt_imp]
gha.dt[is.na(L1_speihishock), L1_speihishock := L1_speihishock_imp]
gha.dt[is.na(L1_speiloshock), L1_speiloshock := L1_speiloshock_imp]
gha <- SpatialPolygonsDataFrame(gha, data.frame(gha.dt), match.ID="rn")
# Moran plots of explanatory vars (e.g. Ghana)
for(i in var) {
  moran.plot(gha@data[[i]], w.gha, zero.policy=T,
    xlab=dt2.lbl[i, varLabel], 
    ylab=paste("Spatially Lagged", dt2.lbl[i, varLabel]),
    labels=gha$svyL2Nm
  )
}


# Compare models
m1 <- lm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt)
m1s <- sacsarlm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt, w.gha, zero.policy=T)

# Results
summary(m1)
summary(m1s)

# Label 20 districts at random
rnd <- sample(1:170, 20)
# Compare models
m1 <- lm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt)
m1s <- sacsarlm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt, w.gha, zero.policy=T)
# Results
summary(m1)

Call:
lm(formula = pcexp_ppp_m ~ spei_lt + L1_speihishock + L1_speiloshock, 
    data = gha.dt)

Residuals:
    Min      1Q  Median      3Q     Max 
-118.41  -38.20   -9.39   24.87  288.70 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     165.165      5.291  31.217  < 2e-16 ***
spei_lt         246.757     53.750   4.591 8.92e-06 ***
L1_speihishock  -80.199     37.767  -2.124   0.0353 *  
L1_speiloshock   -8.735     37.650  -0.232   0.8168    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 64.21 on 159 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1306,    Adjusted R-squared:  0.1142 
F-statistic: 7.964 on 3 and 159 DF,  p-value: 5.573e-05
summary(m1s)

Call:
sacsarlm(formula = pcexp_ppp_m ~ spei_lt + L1_speihishock + L1_speiloshock, 
    data = gha.dt, listw = w.gha, zero.policy = T)

Residuals:
      Min        1Q    Median        3Q       Max 
-116.6812  -37.7663   -9.7939   24.8756  288.1350 

Type: sac 
Coefficients: (asymptotic standard errors) 
               Estimate Std. Error z value  Pr(>|z|)
(Intercept)    138.1537    55.0222  2.5109   0.01204
spei_lt        235.3452    57.4579  4.0960 4.204e-05
L1_speihishock -81.2082    36.3581 -2.2336   0.02551
L1_speiloshock  -7.2666    37.0376 -0.1962   0.84446

Rho: 0.16848
Asymptotic standard error: 0.34074
    z-value: 0.49446, p-value: 0.62098
Lambda: -0.16499
Asymptotic standard error: 0.39596
    z-value: -0.41669, p-value: 0.67691

LR test value: 0.14552, p-value: 0.92982

Log likelihood: -907.6312 for sac model
ML residual variance (sigma squared): 3972.9, (sigma: 63.031)
Number of observations: 163 
Number of parameters estimated: 7 
AIC: 1829.3, (AIC for lm: 1825.4)
# Label 20 districts at random
rnd <- sample(1:170, 20)
# simple LM
plot(gha.dt[!is.na(pcexp_ppp_m), pcexp_ppp_m], m1$fitted.values,
  main="Simple LM", 
  xlab="pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock",
  cex=.5, pch=16)
text(gha.dt[!is.na(pcexp_ppp_m)][rnd, pcexp_ppp_m], m1$fitted.values[rnd], 
  labels=gha.dt[!is.na(pcexp_ppp_m)][rnd, svyL2Nm],
  cex=.6, pos=4)

# SARAR
plot(m1s$y, m1s$fitted.values,
  main="SAC/SARAR LM", 
  xlab="pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock",
  cex=.5, pch=16)
text(m1s$y[rnd], m1s$fitted.values[rnd], 
  labels=gha.dt[!is.na(pcexp_ppp_m)][rnd, svyL2Nm],
  cex=.6, pos=4)

LS0tCnRpdGxlOiAiU1NBIFBvdmVydHkgLSBTcGF0aWFsIFJlZ3Jlc3Npb24iCmF1dGhvcjogIklGUFJJL0hhcnZlc3RDaG9pY2UiCmRhdGU6ICJEUkFGVCA3LzMxLzIwMTYiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazogCiAgICBmaWdfaGVpZ2h0OiA2Ci0tLQoKYGBge3IsIHdhcm5pbmc9RiwgY2FjaGU9VH0KCmxpYnJhcnkocmdkYWwpCmxpYnJhcnkoZm9yZWlnbikKbGlicmFyeShkYXRhLnRhYmxlKQpsaWJyYXJ5KHJoYW5kc29udGFibGUpCmxpYnJhcnkoc3BkZXApCmxpYnJhcnkodG1hcCkKCiMgTG9hZCBkaXN0cmljdHMgYW5kIGF0dHJpYnV0ZXMKZzIgPC0gcmVhZE9HUigiLi9kYXRhIiwgInN2eU1hcHNfMjAxNi4wNi4yMl9zYXJhIikKZHQyIDwtIHJlYWQuZHRhKCIuL2RhdGEvU1NBcG92ZXJ0eV9EaXN0X2ZvckdXUi4xMi5kdGEiKQoKIyBLZWVwIFNUQVRBIGxhYmVscwpkdDIubGJsIDwtIGRhdGEudGFibGUodmFyQ29kZT1uYW1lcyhkdDIpLCB2YXJMYWJlbD1hdHRyKGR0MiwgInZhci5sYWJlbHMiKSkKc2V0a2V5KGR0Mi5sYmwsIHZhckNvZGUpCmR0MiA8LSBkYXRhLnRhYmxlKGR0MikKCiMgTWVyZ2UgYXR0cmlidXRlcyB0byBzaGFwZWZpbGUKZzIuZHQgPC0gZGF0YS50YWJsZShnMkBkYXRhKQpnMi5kdFssIHJuMjIgOj0gYXMuY2hhcmFjdGVyKHJuKV0KZzIuZHRbLCBybiA6PSByb3cubmFtZXMoZzIpXQoKIyBSZWNvZGUgRXRoaW9waWEgd29yZWRhcwpkdDJbc3Z5Q29kZT09ImV0aDIwMTAiLCBzdnlMMkNkIDo9IHN2eUwxQ2QgKiAxMDAwMCArIHN2eUwyQ2RdCgpzZXRrZXkoZzIuZHQsIElTTzMsIHN2eUNvZGUsIHN2eUwxQ2QsIHN2eUwyQ2QpCnNldGtleShkdDIsIElTTzMsIHN2eUNvZGUsIHN2eUwxQ2QsIHN2eUwyQ2QpCgojIEFueSB1bm1hdGNoZWQ/CmR0MlshZzIuZHRdWywgLk4sIGJ5PXN2eUNvZGVdCgojIFJlbW92ZSBgcm5gIGZyb20gU2FyYSdzIHZhcnMKZHQyWywgcm4gOj0gTlVMTF0KCiMgU2VlbXMgb2theSwgc28gbWVyZ2UKZzIuZHQgPC0gZHQyW2cyLmR0XQoKIyBDaGVjayBuYW1lcwpuYW1lcyhnMi5kdClbbmFtZXMoZzIuZHQpICVsaWtlJSAicm4iXQoKCgpgYGAKCmBgYHtyIHNwYXQsIHdhcm5pbmc9Rn0KCiMgTmVpZ2hib3IgbGlzdApuYjIgPC0gcG9seTJuYihnMiwgcm93Lm5hbWVzPXBhc3RlKGcyJElTTzMsIGcyJHJuLCBzZXA9Ii4iKSkKbmIyCgpgYGAKCgpQbG90IGNvbnRpZ3VpdGllcyBmb3IgZS5nLiBHaGFuYQoKCmBgYHtyfQojIFBsb3QgY29udGlndWl0aWVzIChlLmcuIEdIQSkKZ2hhIDwtIGcyW2cyJElTTzM9PSJHSEEiLF0KZ2hhLmR0IDwtIGcyLmR0W0lTTzM9PSJHSEEiXQpnaGEgPC0gU3BhdGlhbFBvbHlnb25zRGF0YUZyYW1lKGdoYSwgZGF0YS5mcmFtZShnaGEuZHQpLCBtYXRjaC5JRD0icm4iKQpjb29yZHMuZ2hhIDwtIGNvb3JkaW5hdGVzKGdoYSkKbmIyLmdoYSA8LSBwb2x5Mm5iKGdoYSkKCnBsb3QoZ2hhLCBib3JkZXI9ImdyZXkiLAogIG1haW49bGlzdCgiQ29udGlndWl0eSAtIEdoYW5hICgxNzAgZGlzdHJpY3RzKSIpKQpwbG90KG5iMi5naGEsIGNvb3Jkcy5naGEsIGNvbD0iYmx1ZSIsIGFkZD1UKQpsZWdlbmQoImJvdHRvbWxlZnQiLCAKICBsZWdlbmQ9Y2FwdHVyZS5vdXRwdXQoc3VtbWFyeShuYjIuZ2hhKSlbMTo3XSwgY2V4PS42KQoKCmBgYAoKYGBge3IsIGNhY2hlPVR9CgojIFNhdmUgZGlzdGFuY2UgbWF0cml4IChXPXJvdyBzdGFuZGFyZGl6ZWQpIHRvIFNUQVRBCncyIDwtIG5iMm1hdChuYjIsIHN0eWxlPSJXIiwgemVyby5wb2xpY3k9VCkKdzIgPC0gYXMuZGF0YS5mcmFtZSh3MikKYXR0cih3MiwgInZhci5sYWJlbHMiKSA8LSBwYXN0ZShnMiRzdnlDb2RlLCBnMiRzdnlMMUNkLCBnMiRzdnlMMkNkLCBzZXA9Ii4iKQp3cml0ZS5kdGEodzIsICIuL2RhdGEvU1NBcG92ZXJ0eV9jb250aW5ndWl0eS5kdGEiLCB2ZXJzaW9uPTEyKQoKYGBgCgoKYGBge3J9CgojIFNwYXRpYWwgd2VpZ2h0cyBmb3IgU1NBIChhbmQgR2hhbmEgb25seSkKI3cyIDwtIG5iMmxpc3R3KG5iMikKdy5naGEgPC0gbmIybGlzdHcobmIyLmdoYSkKCiMgRXhwbGFuYXRvcnkgdmFycwp2YXIgPC0gYygic3BlaV9sdCIsICJMMV9zcGVpaGlzaG9jayIsICJMMV9zcGVpbG9zaG9jayIpCgojIFBsb3QgdG8gaWRlbnRpZnkgZGlzdHJpY3RzIHdpdGggbWlzc2luZyB2YWx1ZXMgKEdoYW5hIG9ubHkpCnRtYXBfbW9kZSgidmlldyIpCnF0bShnaGEsIGZpbGw9dmFyWzFdKQoKIyBJbXB1dGUgbWlzc2luZyB2YWx1ZXMgd2l0aCByZWdpb25hbCBtZWRpYW4KdG1wIDwtIGdoYS5kdFssIGxhcHBseSguU0QsIG1lZGlhbiwgbmEucm09VCksIGJ5PXN2eUwxQ2QsIC5TRGNvbHM9dmFyXQpzZXRuYW1lcyh0bXAsIDI6NCwgcGFzdGUobmFtZXModG1wKSwgImltcCIsIHNlcD0iXyIpKQpzZXRrZXkodG1wLCBzdnlMMUNkKQpzZXRrZXkoZ2hhLmR0LCBzdnlMMUNkKQpnaGEuZHQgPC0gdG1wW2doYS5kdF0KCmdoYS5kdFtpcy5uYShzcGVpX2x0KSwgc3BlaV9sdCA6PSBzcGVpX2x0X2ltcF0KZ2hhLmR0W2lzLm5hKEwxX3NwZWloaXNob2NrKSwgTDFfc3BlaWhpc2hvY2sgOj0gTDFfc3BlaWhpc2hvY2tfaW1wXQpnaGEuZHRbaXMubmEoTDFfc3BlaWxvc2hvY2spLCBMMV9zcGVpbG9zaG9jayA6PSBMMV9zcGVpbG9zaG9ja19pbXBdCgpgYGAKCgpgYGB7cn0KCiMgTW9yYW4gcGxvdHMgb2YgZXhwbGFuYXRvcnkgdmFycyAoZS5nLiBHaGFuYSkKZ2hhIDwtIFNwYXRpYWxQb2x5Z29uc0RhdGFGcmFtZShnaGEsIGRhdGEuZnJhbWUoZ2hhLmR0KSwgbWF0Y2guSUQ9InJuIikKCmZvcihpIGluIHZhcikgewogIG1vcmFuLnBsb3QoZ2hhQGRhdGFbW2ldXSwgdy5naGEsIHplcm8ucG9saWN5PVQsCiAgICB4bGFiPWR0Mi5sYmxbaSwgdmFyTGFiZWxdLCAKICAgIHlsYWI9cGFzdGUoIlNwYXRpYWxseSBMYWdnZWQiLCBkdDIubGJsW2ksIHZhckxhYmVsXSksCiAgICBsYWJlbHM9Z2hhJHN2eUwyTm0pCn0KCmBgYAoKCgpgYGB7cn0KCiMgQ29tcGFyZSBtb2RlbHMKbTEgPC0gbG0ocGNleHBfcHBwX21+c3BlaV9sdCtMMV9zcGVpaGlzaG9jaytMMV9zcGVpbG9zaG9jaywgCiAgZGF0YT1naGEuZHQpCm0xcyA8LSBzYWNzYXJsbShwY2V4cF9wcHBfbX5zcGVpX2x0K0wxX3NwZWloaXNob2NrK0wxX3NwZWlsb3Nob2NrLCAKICBkYXRhPWdoYS5kdCwgdy5naGEsIHplcm8ucG9saWN5PVQpCgojIFJlc3VsdHMKc3VtbWFyeShtMSkKc3VtbWFyeShtMXMpCgpgYGAKCgoKYGBge3J9CgojIExhYmVsIDIwIGRpc3RyaWN0cyBhdCByYW5kb20Kcm5kIDwtIHNhbXBsZSgxOjE3MCwgMjApCmBgYAoKYGBge3J9CgojIFNpbXBsZSBMTQpwbG90KG0xQG1vZGVsJHBjZXhwX3BwcF9tLCBtMSRmaXR0ZWQudmFsdWVzLAogIG1haW49IlNpbXBsZSBMTSIsIAogIHhsYWI9InBjZXhwX3BwcF9tfnNwZWlfbHQrTDFfc3BlaWhpc2hvY2srTDFfc3BlaWxvc2hvY2siLAogIGNleD0uNSwgcGNoPTE2KQoKdGV4dChtMUBtb2RlbCRwY2V4cF9wcHBfbVtybmRdLCBtMSRmaXR0ZWQudmFsdWVzW3JuZF0sIAogIGxhYmVscz1naGEuZHRbIWlzLm5hKHBjZXhwX3BwcF9tKV1bcm5kLCBzdnlMMk5tXSwKICBjZXg9LjYsIHBvcz00KQoKIyBTQVJBUgpwbG90KG0xcyR5LCBtMXMkZml0dGVkLnZhbHVlcywKICBtYWluPSJTQUMvU0FSQVIgTE0iLCAKICB4bGFiPSJwY2V4cF9wcHBfbX5zcGVpX2x0K0wxX3NwZWloaXNob2NrK0wxX3NwZWlsb3Nob2NrIiwKICBjZXg9LjUsIHBjaD0xNikKCnRleHQobTFzJHlbcm5kXSwgbTFzJGZpdHRlZC52YWx1ZXNbcm5kXSwgCiAgbGFiZWxzPWdoYS5kdFshaXMubmEocGNleHBfcHBwX20pXVtybmQsIHN2eUwyTm1dLAogIGNleD0uNiwgcG9zPTQpCgpgYGAKCg==